最近一篇Science的图表中第一幅主图使用圆堆积图(circle packing chart，是树图Treemap/Maptree的一种变体）展示，让我们惊艳到了。通过宏基因组提取16S rRNA基因序列注释细菌群落并统计不同门类细菌差异，最大的圈代表门水平，逐渐缩小的圈按照梯度分别代表纲，科，属。不同颜色标记差异的微生物。这里带大家一起探索这个图表的奥秘。
其实Maptree的中文介绍似乎不多，基于R语言绘制相关图表资料较少。本次教程我将整个函数封装起来，并且定义了三个可简单使用，实现复杂格式转换，一键完成绘制maptreeR函数，相信大家可以很轻易的绘制出Science水准的maptree。
是不是很漂亮
# rm(list=ls())
#导入功能函数
source("./MaptreeForMicrobiome.R")输入文件
文件夹内要准备至少两个文件：OTU表和物种注释
测试数据和代码下载，请在后台回复“maptree”获取下载链接
OTU表otutab.txt格式如下：行为特征OTU/ASV，列为样本名，可以为原始值或标准化的小数均可
OTUID KO1 KO2 KO3
ASV_1 1113 1968 816
ASV_2 1922 1227 2355
ASV_3 568 460 899
物种注释taxonomy.txt：包括OTUID和7级注释，末知的补Unassigned
OTUID Kingdom Phylum Class Order Family Genus Species
ASV_1 Bacteria Actinobacteria Actinobacteria Actinomycetales Thermomonosporaceae Unassigned Unassigned
ASV_2 Bacteria Proteobacteria Betaproteobacteria Burkholderiales Comamonadaceae Pelomonas Pelomonas_puraquae
ASV_3 Bacteria Proteobacteria Gammaproteobacteria Pseudomonadales Pseudomonadaceae Rhizobacter Rhizobacter_bergeniae开始运行
#导入otu表格
otu = read.delim("./otutab.txt",row.names = 1)
head(otu)
otu = as.matrix(otu)
# str(otu)
#导入注释文件
tax = read.delim("./taxonomy.txt",row.names = 1)
head(tax)
#将物种注释文件转化为矩阵，方便使用phyloseq封装
tax = as.matrix(tax)初步绘制maptree图形
data_to_maptree 函数功能为将提供的文件转换为maptree需要的格式，并提供简单出图，参数N代表选取OTU的数量（按照丰度排序，选取丰度最高的前N个OTU做展示），这里我设置N = 200。
在 data_to_maptree 函数中我将几个结果做了封到了一个list中作为输出，
list(p,deg,vertices_t,ps_sub,mygraph)：
p：
maptree图形
deg：
maptree从属关系矩阵
vertices_t：
maptree中OTU以及每层物种的全部信息矩阵
ps_sub：
phyloseq封装的做maptree的OT u信息和物种注释信息。
mygraph：
这是做maptree的网络文件，我将其作为一个输出可以让我我们很容易得到其他ggraph其他出图风格的图形。
参见补充一
#基本函数，画出maptree
mapdata = data_to_maptree (otu,tax,200)
#提取图形部分并保存
p1= mapdata[[1]]
p1
ggsave("./maptree1.pdf", p1, width = 12, height =10 )
按照物种分类信息上色
maptree_add1_plot函数承接data_to_maptree 函数，并使用data_to_maptree 函数的输出作为输入文件，并对data_to_maptree的结果进行按照门水平上色。
#按照平均丰度修改大小和按照门水平上色颜色
mapadd = maptree_add1_plot(mapdata)
p2 = mapadd[[1]]
p2
#d导出ggplot对象的优点就是可以随心所欲的修改图形上的内容
p2 +scale_fill_brewer()
# p2 + scale_color_gradient2()
ggsave("./maptree.pdf", p2, width = 12, height =10)按照OTU差异进行上色
这里我们导入OTU差异分析完成后的表格
#差异表达文件
addtab = read.delim("./KO-WT_all.txt",row.names = 1)
head(addtab)maptree_add2_plot函数承接data_to_maptree 函数，但是填充颜色已经从门变成了自定义列，这也将大大方便我们的实际研究，例如：使用差异分析结果进行上色。
#按照指定列上色
mapadd2 = maptree_add2_plot(mapdata,fillcolor = "level")
#提取作图文件
p3 = mapadd2[[1]]
p3
# p3 +scale_fill_brewer()
ggsave("./maptree3.pdf", p3, width = 12, height =10)
补充一 ：ggraph版本的物种分类树
我们造的轮子在ggraph包中还可以做一些有意思的图形，这里我就本次数据结合ggraph的layout做一个简单的优化展示，为大家带来ggraph版本的物种分类树
graph = mapadd[[2]]
#展示物种分类关系
ggraph(graph, 'dendrogram', circular = TRUE) +
geom_edge_elbow() +
scale_edge_colour_distiller(palette = "RdPu") +
geom_node_point(aes( x = x, y=y, filter = leaf,size=mean,colour=Phylum, alpha=0.2)) +
geom_node_text(aes(x = x*1.05, y=y*1.05, filter = leaf, label=name,
hjust='outward', angle = -((-node_angle(x, y)+90)%%180)+90, size=3,
colour=Phylum), size=0.8, alpha=1) + theme_void()
graph = mapadd[[2]]
# data = create_layout(graph, layout = 'dendrogram', circular = TRUE)
# head(data)
#展示为无数种分类树，这种方式也是我知道
ggraph(graph, layout = 'dendrogram', circular = TRUE) +
geom_edge_diagonal(colour="blue") +
scale_edge_colour_distiller(palette = "RdPu") +
geom_node_point(aes( x = x, y=y, filter = leaf,size=mean,colour=Phylum, alpha=0.2)) +
geom_node_text(aes(x = x*1.05, y=y*1.05, filter = leaf, label=name,
hjust='outward', angle = -((-node_angle(x, y)+90)%%180)+90, size=3,
colour=Phylum), size=0.5, alpha=1) + theme_void()
这是我运行的环境
sessionInfo()
## R version 3.6.1 (2019-07-05)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 18362)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=Chinese (Simplified)_China.936
## [2] LC_CTYPE=Chinese (Simplified)_China.936
## [3] LC_MONETARY=Chinese (Simplified)_China.936
## [4] LC_NUMERIC=C
## [5] LC_TIME=Chinese (Simplified)_China.936
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] RColorBrewer_1.1-2 ggtree_1.14.6 phyloseq_1.28.0
## [4] data.tree_0.7.11 viridis_0.5.1 viridisLite_0.3.0
## [7] forcats_0.4.0 stringr_1.4.0 dplyr_0.8.3
## [10] purrr_0.3.2 readr_1.3.1 tidyr_0.8.3
## [13] tibble_2.1.3 tidyverse_1.2.1 ggraph_2.0.0
## [16] ggplot2_3.2.1.9000 igraph_1.2.4.1
##
## loaded via a namespace (and not attached):
## [1] nlme_3.1-140 lubridate_1.7.4 httr_1.4.0
## [4] tools_3.6.1 backports_1.1.5 vegan_2.5-5
## [7] R6_2.4.1 lazyeval_0.2.2 mgcv_1.8-28
## [10] BiocGenerics_0.30.0 colorspace_1.4-1 permute_0.9-5
## [13] ade4_1.7-13 withr_2.1.2 tidyselect_0.2.5
## [16] gridExtra_2.3 compiler_3.6.1 cli_1.1.0
## [19] rvest_0.3.5 Biobase_2.44.0 xml2_1.2.0
## [22] labeling_0.3 scales_1.1.0 digest_0.6.23
## [25] rmarkdown_1.15 XVector_0.24.0 pkgconfig_2.0.2
## [28] htmltools_0.3.6 rlang_0.4.2 readxl_1.3.1
## [31] rstudioapi_0.10 farver_2.0.1 generics_0.0.2
## [34] jsonlite_1.6 magrittr_1.5 biomformat_1.12.0
## [37] Matrix_1.2-17 Rcpp_1.0.3 munsell_0.5.0
## [40] S4Vectors_0.22.0 Rhdf5lib_1.6.0 ape_5.3
## [43] lifecycle_0.1.0 stringi_1.4.3 yaml_2.2.0
## [46] MASS_7.3-51.4 zlibbioc_1.30.0 rhdf5_2.28.0
## [49] plyr_1.8.4 grid_3.6.1 parallel_3.6.1
## [52] ggrepel_0.8.1 crayon_1.3.4 lattice_0.20-38
## [55] splines_3.6.1 graphlayouts_0.5.0 Biostrings_2.52.0
## [58] haven_2.1.1 multtest_2.40.0 hms_0.5.0
## [61] zeallot_0.1.0 knitr_1.23 pillar_1.4.2
## [64] reshape2_1.4.3 codetools_0.2-16 stats4_3.6.1
## [67] glue_1.3.1 evaluate_0.14 data.table_1.12.2
## [70] modelr_0.1.4 treeio_1.6.2 vctrs_0.2.0
## [73] tweenr_1.0.1 foreach_1.4.4 cellranger_1.1.0
## [76] gtable_0.3.0 polyclip_1.10-0 assertthat_0.2.1
## [79] xfun_0.9 ggforce_0.3.1 broom_0.5.2
## [82] tidygraph_1.1.2 tidytree_0.2.4 survival_2.44-1.1
## [85] rvcheck_0.1.3 iterators_1.0.10 IRanges_2.18.1
## [88] cluster_2.1.0撰文：文涛，南京农业大学
责编：刘永鑫，中科院遗传发育所
